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Abstract 

Direct numerical simulations (DNS) of flow over an NACA-0012 airfoil are performed at a low and a 
moderate Reynolds numbers of Re,=50*10° and 1*10°. The angles of attack are 5 and 15 degrees at the 
low and the moderate Reynolds number cases respectively. The three-dimensional unsteady 
compressible Navier-Stokes equations are solved using higher order compact schemes. The flow field in 
the low Reynolds number case consists of a long separation bubble near the leading-edge region and an 
attached boundary layer on the aft part of the airfoil. The shear layer that formed in the separated 
region persisted up to the end of the airfoil. The roles of the turbulent diffusion, advection, and 
dissipation terms in the turbulent kinetic-energy balance equation change as the boundary layer 
evolves over the airfoil. In the higher Reynolds number case, the leading-edge separation bubble is 
very small in length and in height. A fully developed turbulent boundary layer is observed in a short 
distance downstream of the reattachment point. The boundary layer velocity near the wall gradually 
decreases along the airfoil. Eventually, the boundary layer separates near the trailing edge. The 
Reynolds stresses peak in the outer part of the boundary layer and the maximum amplitude also 
gradually increases along the chord. 


I. Introduction 


Flow over a two-dimensional airfoil exhibits different flow dynamics and aerodynamic characteristics depending 
on the shape of the airfoil and the freestream parameters.” The shape of an airfoil is defined by the camber, 
thickness distribution, leading edge radius, and trailing edge thickness. The freestream parameters that influence the 
flow dynamics are the Mach number, Reynolds number, and angle of attack. The shape and the angle of attack 
basically determine the pressure distribution on the airfoil. The pressure distribution on the pressure side or the 
lower surface of the airfoil is favorable, and the pressure distribution on the suction side or the upper surface is 
favorable near the leading-edge region and becomes adverse on the downstream part of the airfoil. The adverse 
pressure gradient induces early boundary layer transition and/or flow separation near the trailing edge part of the 
airfoil. The leading-edge radius and angle of attack also can promote flow separation near the leading edge of the 
airfoil. These separations, near the leading edge and the trailing edge, are strongly influenced by the Reynolds 
number and the angle of attack. 

At low Reynolds numbers and at small angles of attack, the flow remains laminar near the leading edge and 
separates due to the strong adverse pressure gradient that exists near the leading-edge region. The separated region 
grows downstream and forms a strong shear layer. This shear layer becomes unstable and breaks down into a 
turbulent state. The turbulent separated shear layer reattaches on the airfoil and forms a long separation bubble on 
the surface of the airfoil. With increasing angles of attack, the shear layer does not reattach on the airfoil and a large 
separated flow forms over the entire upper surface. At higher Reynolds numbers, the leading-edge separation bubble 
becomes small and an attached turbulent flow prevails downstream up to the trailing edge. With increasing angle of 
attack, the flow starts to separate near the trailing edge. The trailing edge separation moves upstream and the 
separated flow grows in size with increasing angle of attack. At one point, the leading edge and the trailing edge 
separated region merge and the flow becomes completely separated over the upper surface of the airfoil. 

The required computer resources and time constraints hinder the solution of the Navier-Stokes (N-S) equations 
numerically for turbulent flows at high Reynolds numbers.’ Hence, the current state of the art is to solve some 
approximate versions of the N-S equations that require modeling. The simplest and the most popular approximate 
method is the Reynolds-Averaged Navier-Stokes (RANS) equations, where the equations are derived for time- 
averaged quantities.* The unclosed Reynolds stresses are modeled to close the equations.’ In general, existing 
models provide good results compared to experiments in attached flows where turbulence is in equilibrium. 
However, prediction of separated flows with existing RANS turbulence models is still not very satisfactory. ’ It is 


difficult to accurately capture both the separation and reattachment points (i.e., the size of the separation bubble). 
The second approximate method uses the Large-Eddy-Simulation (LES) equations, which are obtained for the large- 
scale motions by filtering small scales. The unknown stresses introduced by the small scales (subgrid-scale stresses, 
or SGS) are modeled to close the equations.® This approach is computationally more expensive than the RANS 
approach, but yields more accurate results without many ad-hoc modifications to the SGS models. This 
methodology is being applied to a wide range of problems. However, these simulations also become very expensive 
at flight Reynolds numbers due to the required grid resolutions to resolve the turbulent structures near the wall 
region in wall bounded turbulent flows. A feasible approach to compute turbulent flows at high Reynolds numbers 
is to employ hybrid methods where the wall region is modeled using some form of RANS equations and the outer 
regions are simulated using LES. 

The first successful direct numerical simulation (DNS) was performed for a channel flow at a Reynolds number, 
based on the bulk velocity and the channel half-width, of 3000.” With increasing computational capability, the 
application of DNS to more complex problems at high Reynolds numbers is currently being pursued by many 
researchers. The early status of DNS in turbulent flow is reviewed by Moin and Mahesh.'° Recently, DNS for a 
channel flow has been extended to a Reynolds number of 125000.'' DNS of a turbulent boundary layer over a flat 
plate was first performed up to a Reynolds number based on momentum thickness of Reg= 1410.'* The computed 
turbulent statistical quantities agreed fairly well with the experiments and the available theories. The simulations 
over a flat plate at zero pressure gradient have been recently extended to higher Reynolds numbers. Schlatter et al.’ 
has performed a DNS simulation over a flat plate up to a Reynolds number of Reg = 2500. The results were 
compared with the detailed experimental results obtained at high Reynolds numbers.'* The finding was that at high 
Reynolds numbers, the wall layer dynamics scale with inner viscous scales and with the outer large-scales when 
compared to the observations at low Reynolds numbers. Several simulations and detailed experiments have been 
performed over a flat plate in adverse pressure gradients.'>? It was found that large-scale motions in the outer part 
of the boundary layer are more energized with increasing pressure gradients. The turbulent intensities first peak in 
the near wall region and also show a broad distribution with another peak in the outer region. 

Breakdown of laminar separation bubbles to turbulence has been simulated numerically in several 
investigations.°”” The separation is forced on a flat plate boundary layer by applying suction through an upper wall. 
The separated shear layer transitions to turbulence, due to the Kelvin-Helmholtz instability, and reattaches and 
evolves downstream as a turbulent boundary layer on a flat plate. The turbulent characteristics in the separation 
region and in the vicinity of the reattachment region resemble the characteristics of a plane mixing layer. It was also 
observed that the evolving boundary layer approaches the equilibrium boundary layer over a flat plate at several 
bubble lengths downstream of reattachment. DNS of flow over airfoils have been performed at low Reynolds 
numbers at low and high angles of attack.*”> At low angles of attack, the flow separates near the leading edge and 
reattaches as a turbulent flow downstream. With increasing angles of attack, the flow does not reattach and remains 
separated over the airfoil. Not many DNS have been performed at high Reynolds numbers due to the required 
computational resources. Typically, LES are performed to predict the flow features over airfoils at high angles of 
attack. In a Large-Eddy-Simulation of flow around the Aerospatiale A-profile*” at a high angle of attack near stall, 
it was found, comparing with experiments, that the treatment of the laminar separation bubble is essential in 
predicting the aerodynamic quantities and the characteristics of the boundary layer in the aft part of the airfoil. 

The objective of the present research is to perform Direct Numerical Simulations (DNS) of flow over an NACA- 
0012 airfoil at low and moderate Reynolds numbers and at different angles of attack to investigate the flow 
dynamics, the turbulent statistical quantities, and the turbulent kinetic energy balance. The data will be used to 
validate and improve future RANS and LES for these cases. The Reynolds numbers based on the freestream 
velocity, density, and the chord length simulated in this paper are Rez = Pudol / Uo = 50*10° and 1*10°. The 
computations were performed at 5 degrees angle of attack at the low Reynolds number case and at 15 degrees for the 
higher Reynolds number case. The freestream Mach number for the low and higher Reynolds numbers cases are M = 
0.4, and 0.2, respectively. The low Reynolds number case was simulated by Jones et al. (2008) and those results 
will be compared with the results of the current study. The computations are performed using 6th-order accurate 
compact scheme for spatial discretization and a 3rd-order Runge-Kutta scheme for time discretization. 


II. Models and Flow Conditions 


Computations are performed for a flow over an NACA-0012 airfoil. The NACA-0012 airfoil with a sharp 
trailing edge is defined by the following equation”® 


y(x) = Cy (coVx + c3x + C4x? + €5x3 + c6x*) (1) 


where 
Cc, = 0.594689180 
C2 = 0.298222773 
c3 = —0.127125232 
C4 = —0.357907906 
Cs = 0.291984971 
C6 = —0.105174606 (2) 


The chord length, c, is used to nondimensionalize the coordinates. The airfoil defined using this equation 
has a zero thickness at the trailing edge. The compact scheme that is used in this work needs a smooth 
grid around the airfoil. The sharp trailing edge was replaced with a circular trailing edge having a small 
radius of 0.001. The modified and the unmodified airfoils are shown in Fig. 1. The new chord length of 
the airfoil is 0.983. However, the length scales are continued to be nondimensionalized by the chord 
length of the unmodified airfoil in this paper. A body-fitted O-type curvilinear grid system was used in all 
of the simulations. The Cartesean coordinate system with the x-axis along the chord, y-axis perpendicular 
to the chord, and the z-axis along the span was also used. The corresponding velocity components are 
denoted by u, v, and w, respectively. The variables density p, pressure p, and the velocities are non- 
dimensionalized by their corresponding freestream variables p.., p, and U. ,respectively. The time 
is nondimensionalized by (c/U..). The computational domain and the grid system are depicted in Fig. 2. 
The leading edge of the airfoil is located at (x, y) = (0, 0). The outer boundary follows a circle of radius R, 
centered at the midchord of the airfoil. The grid is uniform in the spanwise z-direction. 
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Figure 1. NACA-0012 Airfoil with a sharp and a rounded trailing edges. 
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Figure 2. Flow over a NACA-0012 Airfoil (a) computational domain, (b) grid distribution (every 10° points are shown). 


Simulations are performed at two chord Reynolds numbers and at different angles of attack. Table 1 gives the 
Reynolds numbers, angles of attack, and freestream Mach numbers considered in the simulations. The Reynolds 
numbers based on the freestream velocity, density, and the chord, simulated in this paper are Reg = PoUVac / Ho = 
50*10°, and 1*10°. The computations are performed at 5 degrees angle of attack for the low Reynolds number case 
and at 15 degrees for the higher Reynolds number case. The freestream Mach number for the low and higher 
Reynolds numbers cases are M = 0.4, and 0.2, respectively. 


Table 1 Freestream parameters used in the simulations. 


Re, Angle of attack, a (deg) Mach number 
50*10° 5 0.4 
1*10° 15 0.2 


ll. Governing Equations 


The partial differential equations solved are the three-dimensional unsteady compressible Navier-Stokes 
equations in conservation form. ’ The viscosity is computed using Sutherland’s law and the thermal conductivity is 
obtained from the Prandtl number. A constant value of 0.70 is assumed for the Prandtl number. A body-fitted 
curvilinear grid system was used in all of the simulations. The equations are transformed from the physical 
coordinate system (x, y, z) to the computational curvilinear coordinate system (€, n, C). 


A. Solution Algorithm 

The governing equations were solved using higher order compact schemes with domain decomposition. 
Flow simulations were performed using a 6th-order compact scheme in all three directions. The orders are reduced 
to 4" and 3" orders near the walls. A higher order implicit filtering” of 8" order was used in all three directions to 
remove the high wavenumber contents in the solution. A 3"order total-variation-diminishing (TVD) Runge-Kutta 
scheme was employed for time integration. The no-slip and constant temperature conditions are imposed at the 
airfoil surface. The wall temperature is fixed at a constant temperature of 300K. Periodic conditions are applied in 
the spanwise direction. Characteristic boundary conditions are used at the outer boundary. A spanwise domain size 
of 0.2c and an outer boundary of radius R, = 15c were selected for the simulations. 
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IV. Results 


A. DNS at Reynolds number Re, = 50*10° and a = 5° 
The DNS results at the lower Reynolds number of Re.= 50*10° at an angle of attack of 5 degrees are presented 
first in this section. 


(1) Two-dimensional simulations 

The two-dimensional simulations were performed with different outer boundary domains R, and grid sizes 
around the airfoil to determine the effect of the outer boundary on the pressure distributions on the airfoil. The 
simulations were conducted with the domain sizes of R. = 15, 30, and 100. The number of grid points around the 
airfoil were varied from 1501 to 3001. There are 501 grid points in the normal direction. Figures 3(a) and (b) depict 
the mean pressure and skin friction coefficients obtained with different domain sizes and grid distributions around 
the airfoil. It is seen that the results obtained with large and small domain sizes are almost the same except small 
differences in the pressure coefficients near the reattachment region. The differences between the finer and coarser 
grids are also negligible. Hence in the subsequent three-dimensional simulations, a domain size of R. = 15 and 
(3001*501) points around the airfoil and in the normal direction were used. The boundary layer on the upper surface 
separates near the leading edge of the airfoil at x/c = 0.135 and reattaches near the midchord at x/c = 0.555. These 
locations are close to those reported in Ref. 23 (i.¢., x/c = 0.151 and 0.582, respectively). 

Figures 4(a) and (b) show the contours of the spanwise vorticity at a fixed time. Figure 4(a) depicts the contours 
over the airfoil and in the wake region and Fig. 4(b) depicts the contours near the reattachment region of the airfoil. 
The separated shear layer rolls up into clockwise isolated vortices between x/c ~ 0.40 and 0.47. Near x/c ~ 0.5 a 
counterclockwise vorticity region protrudes from the surface into the evolving shear layer. Downstream of this 
region, x/c ~ 0.5, isolated clockwise rotating vortices are shed. Figure 5(a) shows the pressure perturbations at a 
fixed time along the upper and lower surfaces of the airfoil. It is to be noted that no external disturbances were 
introduced in the simulations. All the unsteadiness is originating naturally from some numerical imbalances and/or 
from some self-sustaining mechanisms. From the airfoil leading edge up to x/c ~ 0.4, the pressure perturbations 
remain small on the order of ~ 10°. The perturbations grew to about 10% of the free stream value within a short 
distance from x/c ~ 0.45 to 0.55. The maximum pressure perturbations after reattachment are about 5% of the 
freestream value. Figure 5(b) shows the variations of total lift (C.) and drag coefficients (Cp,r) with time. The drag 
contributions from pressure (Cp) and skin friction (Cpr) are also included separately in the figure. It is seen that lift 
and drag coefficients are not exactly periodic in time. The approximate period of oscillations is about 0.27 and the 


frequency of shedding is about 3.7. The time averaged lift and drag are 0.502 and 0.0315, respectively. The 
contributions from the frictional and pressure drag are 0.0095 and 0.022, respectively. Even at this small angle of 


attack of 5 degrees, the pressure-induced drag contributes to about 70% of the total drag and the skin friction 
contributes to only about 30%. 
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Figure 4. Contours of spanwise vorticity at a fixed time (a) over the airfoil and in the wake and (b) near the 
reattachment region of the airfoil. 


(a) (b) 


Lower 
Upper 


> 0.05 + | 

2 

Qa 

~ 0.00 F i) | 
0.05 F | 


0.10 + 


0.0 0.2 0.4 0.6 0.8 7.0 0.0 05 7.0 75 2.0 25 3.0 
xlc t 


Figure 5. (a) Pressure perturbations over the airfoil at a fixed time and (b) variations of lift and drag with time. 


(2) Three-dimensional simulations 

The three-dimensional direct numerical simulations (DNS) were performed with a spanwise domain size of 0.2c. 
This was the domain size used in the simulations by Jones et al.”? The simulations were performed with grid sizes of 
(1501*501*145) and (3001*501*201). The numbers in the brackets reveal the number of points in the streamwise, 


normal, and spanwise directions, respectively. They are referred to as coarse and fine grids in the discussions below. 
The grid spacings in viscous units change along the streamwise direction due to the change in friction velocity along 
the airfoil surface. Figure 6 shows the variations of the grid spacing in wall units in the x, y, and z directions along 
the upper surface of the airfoil for the coarse and fine grids. The maximum grid spacing in the x- and z-directions are 
both below 3 for the fine grid and below 6 and 4 for the coarse grid. The minimum grid spacing at the wall varies 
between Ay’ = 0.1 to 0.3. These resolutions are much finer than that used by Jones et al.”? Figures 7(a-c) show the 
variation of the ratio of grid spacing in each direction to the Kolmogorov length scale, 7, along the normal direction 
at different stations, x/c = 0.02, 0.2, 0.3, 0.4, 0.5, 0.6, and 0.7 along the upper surface of the airfoil. The results for 
the fine grid and the coarse grids are displayed at one station x/c = 0.7. The fine grid sizes to Kolmogorov scale ratio 
are below 3 in the streamwise and spanwise directions. This ratio is below 2 for the grid distribution in the normal 
direction. It will be observed later that these grid distributions are sufficient to resolve the flow field up to the 
dissipation scales. 
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Figure 7. Ratio of grid spacing to Kolmogorov scale in the normal direction at different stations along the upper surface 
of the airfoil. Re, = 50*10°, a=5°. 


The transport equations for the Reynolds stresses and the turbulent kinetic energy are given in Ref. 30. 
The equation for the turbulent kinetic energy can be written as 


The left hand side of Eq. (3) is the local rate of change of turbulent kinetic energy and the terms on the right 
hand side are: (1) production, (2) diffusion, (3) dissipation, and (4) advection. In the statistical steady state 


conditions, the local rate of change should be zero and these four terms should balance each other to yield zero sum. 
A nonzero left hand side points to a lack of grid resolution in the simulations. Figures 8(a) and (b) depict the 
magnitudes of the four terms on the right hand side of this equation and the balance or the total of these four terms at 
the stations x/c = 0.6 and 0.8 for the coarse and fine grids. It is seen that the balance is below 10% of the production 
at x/c = 0.6 and 0.8 stations for the fine grid and they are below 15% for the coarse grid. This suggests that the grid 
used in this DNS resolved most of the scales including the dissipation scales at this Reynolds number. The results 
obtained with the fine grid will be presented in the following sections. 
(a) x/c = 0.6 (b) x/c = 0.8 
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Figure 8. Balance of different terms in the turbulent kinetic equation using different grids. (a) x/c = 0.6, (b) x/c = 0.8. Re. = 
50*10°, a=5°. 


Figures 9(a) and (b) display the instantaneous isosurfaces of the second invariant of the velocity gradient tensor 
at O = 500 and the spanwise vorticity contours in the (x, y) plane at the spanwise location z = 0.1, respectively. The 
flow remains steady and laminar up to x ~ 0.4. Close to x ~ 0.45 large scale vortices are shed from the shear layer. 
These vortices break down into turbulence downstream. 


Figure 9. Instantaneous (a) iso-surfaces of Q@ = 500 and (b) spanwise vorticity contours in the (x-y) plane at z = 0.1. Re. = 
50*10°, a=5°. 


The statistical quantities are obtained by averaging the solution in the spanwise direction and in time. The 
averaging in time was performed for about 5 flow periods. Figure 10 shows the mean streamwise velocity contours 
and the streamlines from the current DNS. Figure 10(a) displays the contours and the streamlines over the entire 
airfoil and Fig. 10(b) displays the contours near the separation bubble. A long separation bubble near the leading 
edge of the airfoil and a thick attached shear layer in the trailing edge part of the airfoil are observed. An interesting 
observation is that the bubble size increases gradually up to x ~ 0.5 and reattaches steeply close to x ~ 0.6. Figures 
11(a) and 11(b) show the pressure coefficient c, and the skin friction cyalong the lower and the upper surfaces of the 
airfoil. The results obtained with the coarse grid are included in these figures. There is a small difference in the skin 
friction distribution near the reattachment region, but the overall agreement is very good between the coarse and fine 


grids. The pressure distribution depicts a plateau region from x ~ 0.1 to 0.5, rises steeply until x ~ 0.6, and recovers 
slowly to the trailing edge value. The skin friction curve shows that the flow separates at about x = 0.1 and 
reattaches at x = 0.6. These values are in good agreement with the computed values of 0.09 and 0.6 given in Ref. 23. 
Table 2 shows the drag and the lift coefficients obtained from the current DNS and the DNS results of Ref. 10. The 
individual contribution from the pressure and the friction to the total drag and lift were also computed. The 
computed total lift coefficient of 0.627 agrees very well with the value of 0.621 obtained by Ref. 23. All the lift is 
coming from the pressure and the contribution from the friction is negligible. The computed drag coefficient of 
0.0375 slightly differs from the value of 0.0358 reported in Ref. 23. The form drag is about 76% of the total drag 
and the frictional drag is the remaining 24%. The results show that for flow over airfoils at low Reynolds numbers, 
even at small angles of attack, most of the drag force is coming from the pressure. 
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Figure 10. Mean u-velocity contours and the streamlines. (a) Flow over the airfoil, and (b) flow near the separation 
bubble. Re. = 50*10°, a=5°. 
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Figure 11. Mean (a) pressure, and (b) skin friction coefficients. Re. = 50*10°, a=5°. 


Table 2 Mean drag and Lift coefficients. 


Ca C 
Pressure __ Friction Total Pressure Friction Total 
Present 0.0270 0.0090 0.0354 0.625 0.001 0.626 
Ref. 23 0.0278 0.0081 0.0358 0.621 


As observed earlier, the flow field over the upper surface of the airfoil consists of two parts. One is the separated 
region near the leading edge and the second is the attached flow region downstream of the separated region. The 
boundary layers in the separation bubble consist of reverse flows near the wall and wake-like profiles with 
inflections in the outer part. Downstream of the reattachment point, a new boundary layer develops along the surface 
of the airfoil and the inflectional profile in the outer part becomes weaker. If the evolution distance is long enough, 
the boundary layer profiles will be similar to the profiles over flat surfaces in adverse pressure gradients. 


In Fig. 12, the computed mean velocity <U,> profiles are shown at different locations along the upper surface of 
the airfoil. Here, U; is the velocity parallel to the surface of the airfoil. Figure 12(a) shows the mean velocity profiles 
near the leading edge and Fig. 12(b) shows the profiles farther downstream. It is seen that the boundary layer 
thickness increases steeply in the separation bubble region from 0.004 at x = 0.02 to 0.02 at x = 0.2. The boundary 
layer profile at x = 0.2 consists of a reverse flow near the surface and a strong inflectional profile in the outer part. 
After the flow attaches at x = 0.6, a two-layer structure develops inside the boundary layer. One is a layer very close 
to the wall and the other is the very thick wake-like profile in the outer part of the boundary layer. The outer layer is 
reminiscent of the inflectional profile that formed in the separation bubble. Even at x ~ 0.8, a weak inflectional 
profile was observed in the outer part of the boundary layer. 
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Figure 12. Mean boundary layer velocity profiles, U,, at different stations. 


In Figs. 13-14, the computed Favre-averaged Reynolds stresses <pu”u’> and <pu’v"> profiles are shown at 
different streamwise locations along the upper surface of the airfoil. The results for the Reynolds stresses show that 
fluctuations are very small ( less than ~ 0.05) near the leading edge up to x ~ 0.4. Turbulent intensity and the shear 
stress peak near the reattachment point of x ~ 0.6. These peaks occur in the middle (y direction) of the wake part of 
the boundary layer. The maximum turbulent intensity increased by about 15 times from the middle of the separation 
bubble x ~ 0.4 to the reattachment point x ~ 0.6. The maximum Reynolds stresses continue to decrease from the 
reattachment point with increasing distance along the airfoil. Near the trailing edge of the airfoil, x > 0.8, the normal 
stress has similar amplitudes near the wall and in the middle part of the boundary layer. This implies that away from 
the separation bubble a wall layer similar to flows over flat surfaces is gradually developing over the airfoil surface. 
This behavior is also observed in the Reynolds shear stress profiles in Fig. 14(b). In turbulent boundary layers over a 
flat plate at zero pressure gradient, the Reynolds stresses peak close to the wall inside the buffer region.” 
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Figure 13. Streamwise normal Reynolds stress, <pu"u'>, profiles at different streamwise stations. 
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Figure 14. Reynolds shear stress, <pu”v”>, profiles at different streamwise stations. 


Figures 15(a-d) depict the magnitudes of the four terms on the right hand side of the turbulent kinetic energy 
equation, Eq. (3), and the balance or the total of these four terms at four stations x = 0.4, 0.5, 0.6, and 0.9. The 
balance is below 5% of the production at x = 0.4 and 0.5, and below 10% and 15% of the production at x = 0.6 and 
0.9, respectively. The first observation is that away from the wall turbulent production is the largest at all streamwise 
stations. The turbulent dissipation and diffusion have very large values near the wall. The large dissipation at the 
wall is balanced by the positive turbulent diffusion term. Up to the reattachment point, the dissipation is small 
compared to the advection and turbulent diffusion terms in the outer part of the boundary layer. The production is 
balanced by the advection and the turbulent diffusion terms. Near the reattachment region, x/c = 0.6, the positive 
advection term near the wall region increases and the contribution from the dissipation also increases. Evolving 
downstream toward the trailing edge of the airfoil the dissipation and production have comparable magnitudes. The 
advection and the turbulent diffusion have positive and negative values across the boundary layer, respectively. It is 
also observed that near the trailing edge at x ~ 0.9, the production has a peak near the surface. This is due to the 
developing boundary layer along the surface. Hence the distributions and the roles of advection, turbulent diffusion, 
and dissipation in the balance of turbulent kinetic energy gradually change when the flow evolves from the 
separation region to the attached boundary layer region. In DNS of the breakdown of short laminar bubbles formed 
on a flat plate,” similar observations also were made. It was observed that the flow near the reattachment consists 
of an outer region that resembles a turbulent mixing layer. It was also noted that the developing wall boundary layer 
along the wall recovers slowly toward an equilibrium boundary layer. 
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Figure 15. Balance of different terms in the turbulent kinetic equation. Re, = 50*10°, a=5°. 


The contours of the average statistical quantities- (a) kinetic energy, (b) shear stress, (c) production, and (d) 
dissipation - are plotted in Figs. 16(a-d). These quantities play important roles in RANS modeling and development. 
It is interesting to see that all these terms peak in the middle part of the boundary layer near the reattachment point. 
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Figure 16. Contours of the turbulent statistical quantities (a) Kinetic energy, (b) Shear stress, (c) Production, and (d) 
Dissipation. Re, = 50*10°, a=5°. 
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B. DNS at moderate Reynolds number Re, = 1*10° and a = 15° 

The DNS were performed at a larger Reynolds number of Re, = 1* 10° and at a high angle of attack of a = 15°. 
There were 3001, 501, and 513 points along the airfoil surface, in the wall normal, and spanwise directions, 
respectively. A spanwise domain size of 0.2c and the outer boundary of the computational domain size of R, = 15c 
were selected in for the simulations. Figures 17(a-c) show the variation of the ratio of grid spacing in each direction 
to the Kolmogorov length scale, 7, along the normal direction at different stations, x/c = 0.02, 0.2, 0.3, 0.4, 0.5, 0.6, 
and 0.7 along the upper surface of the airfoil. The grid sizes compared to Kolmogorov scale are below 5 in the 
normal and spanwise directions. This ratio is below 15 for the grid distribution in the streamwise direction. The grid 
distribution in the spanwise direction is sufficient to resolve the flow field up to the dissipation scales. However, the 
grid distributions in the streamwise and normal directions are about two times too coarse to capture the dissipation 
scales accurately. The number of points in the streamwise and normal directions have to be redistributed or 
increased to resolve the dissipation scales. 


0 5 Ayly 10 


Figure 17. Grid spacing in wall units along the upper surface of the airfoil. Re, = 1*10°, a=15°. 


Figure 18(a) shows the instantaneous isosurfaces of the second invariant of the velocity tensor with a value O = 
500. The isosurfaces are colored by the local velocity parallel to the airfoil surface. Figure 18(b) shows the 
instantaneous w-velocity in the plan view (x-z plane) along a fixed normal grid (j=20, y ~ 10), which is located a 
short distance above the lower wall. The white region displays the region with the negative velocity. The 
instantaneous separation regions extend to larger domains compared to those derived from the time averaged flow 
field. This phenomenon is common in separated flows.” The statistical quantities are obtained by averaging the 
solution in the spanwise direction and in time. The averaging in time was performed for about 2 flow periods. Figure 
19 shows the mean streamwise velocity contours and the streamlines from the current DNS. Figure 19(a) displays 
the contours over the entire airfoil and Fig. 19(b) displays the contours near the separation bubble. There exists a 
small separation bubble near the leading edge of the airfoil and a second bubble near the airfoil’s trailing edge. The 
flow separates near x ~ 0.005 and reattaches at x ~ 0.024. The trailing edge separation occurs at x ~ 0.80 and 
reattaches at the trailing edge. No freestream turbulence was included in the simulations. The boundary layer near 
the stagnation region remains laminar. The laminar boundary layer that developed on the upper surface separates 
near the leading edge due to the strong local adverse pressure gradient. Since the Reynolds number is higher, the 
detached shear layer becomes strongly unstable. This unstable shear layer breaks down into turbulence and 
reattaches within a short distance from the separation point. The breakdown process depends on the initial 
disturbances generated in the laminar region of the separation bubble. It was experimentally observed hat high 
freestream turbulence level and artificial tripping near the leading edge strongly influence the lift and drag of 
airfoils. 
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Figure 18. Instantaneous (a) iso-surfaces of Q@ = 500 and are colored by the local u-velocity parallel to the surface (0.2 
<u < 1.2), and (b) contours of u-velocity in the (x-z) plane at y = 0.05. Re, = 1*10°, a=15°. 


(a) Entire airfoil (b) Separation bubble near the leading edge 


Figure 19. Mean u-velocity contours and the streamlines. (a) Flow over the airfoil, and (b) flow near the separation 
bubble. Re, = 1*10°, a=15°. 


Figures 20(a) and (b) display the instantaneous pressure distributions over the airfoil. Figure 20(a) displays the 
pressure distribution over the entire airfoil and Fig. 20(b) displays the distribution near the separation bubble. The 
stagnation point is located on the lower surface at x ~ 0.039. As observed earlier, the pressure distribution on the 
lower surface remains favorable over the airfoil. The flow accelerates from the stagnation point towards the rounded 
leading edge causing the pressure to steeply decrease around the leading edge. The pressure then increases strongly 
from the leading edge causing the flow to separate slightly downstream. The pressure distribution forms a plateau up 
to the middle of the separation bubble and starts to oscillate strongly from x ~ 0.014. At the end of the oscillations 
near x ~ 0.023, the pressure increases steeply. Beyond this point, the pressure increases mildly towards the trailing 
edge. 
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Figure 20. Instantaneous pressure distribution. (a) Over the full airfoil, and (b) near the leading edge. Re, = 1*10°, 
a=15°. 


Figure 21 shows the mean pressure distribution over the airfoil. The RANS calculations were performed using 
CFL3D code* assuming the flow is fully turbulent over the airfoil. The agreement between the DNS and the RANS 
computations are good except over the separation bubble near the leading edge region. In the RANS computation, 
the separation bubble is very narrow compared to the DNS results. This is because in the DNS the flow remains 
laminar up to the middle of the separation bubble, hence the separation bubble is longer in the DNS than in RANS. 
Figures 22(a) and (b) show the skin friction Cy along the lower and the upper surfaces of the airfoil as well as the 
results from the RANS computations. The skin friction on the lower surface differs from the RANS solution because 
the flow remains laminar on the lower surface in the DNS, whereas it is turbulent in the RANS computations. The 
skin friction on the upper surface agrees with the DNS results in the aft part of the airfoil on the upper surface. At 
this angle of attack, the RANS results compare satisfactorily with those of the DNS results in the trailing edge 
region of the airfoil. Table 3 shows the drag and the lift coefficients obtained from the current DNS. The individual 
contribution from the pressure and the friction to the total drag and lift were also calculated. All the lift is coming 
from the pressure and the contribution from the friction is negligible. The computed drag coefficient of 0.0557. The 
form drag is about 92% of the total drag and the frictional drag is the remaining 8%. The results show that for flow 
over airfoils at this high angle of attack, where the flow has a small separation near the trailing edge, most of the 
drag force is coming from the pressure. 


Table 3 Mean drag and Lift coefficients. 


Ca C 
Pressure Friction Total Pressure Friction Total 
0.0515 0.0042 0.0557 1.290 0.000 1.290 


14 


Lower DNS 
Upper 
— -— — - Lower CFL3D sst 
— — — — Upper 


Figure 21. Mean pressure coefficient obtained from the DNS and the RANS computations. Re, = 1*10°, a=15°. 
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Figure 22. Mean skin friction coefficients. Re. = 1*10°, a=15°. 


Figures 23(a) and (b) show the mean boundary layer profiles at different locations along the upper surface of the 
airfoil. Figure 23(a) depicts the profiles inside the separation bubble at x ~ 0.02 and immediately downstream of the 
reattachment point at x ~ 0.04, 0.1 and 0.2. Figure 23(b) displays the results in the attached part of the airfoil at x ~ 
0.4 and 0.6, and in the trailing edge separation bubble region at x ~ 0.8 and 0.9. It is seen that downstream of the 
reattachment point velocity profiles with strong gradients near the wall are developed. When the boundary layer 
evolves downstream, it thickens and the velocity near the wall decreases. It is interesting to note that even before the 
separation point x ~ 0.8, inflectional profiles start to develop at x ~ 0.6. 
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(a) x = 0.02, 0.04, 0.1, and 0.2 (b) x = 0.4, 0.6, 0.8, and 0.9 
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Figure 23. Mean boundary layer velocity profiles at different stations. 

Figures 24(a-f) show all the Favre averaged Reynolds stress components < pu”u" >,< pv'"'v" >, < pw"w" > 
,and < pu"'v" > profiles at six stations x = 0.02, 0.2, 0.4, 0.6, 0.8 and 0.9. The mean velocity profiles were also 
included in these figures. The first station, x = 0.02, is located in the middle of the leading edge separation bubble. 
The stations x = 0.2, 0.4, and 0.6 are in the attached region of the airfoil. The last two stations 0.8, and 0.9 are 
situated inside the separation zone near the trailing edge. The first observation is that the root-mean-square values in 
the leading edge separation bubble are one order of magnitude higher than in the downstream part of the airfoil. In 
the leading edge separated region x = 0.02, the turbulent normal stresses peak in the similar maximum shear region 
as we observed in the low Reynolds number case. It is also observed that shear stress peaks near the surface and the 
normal stresses also have large values near the surface. Immediately downstream of the separation bubble at x = 
0.20, and 0.40, there exist two layers, one very near the surface and the other in the wake part of the boundary layer. 
The thin layer near the wall has large normal stresses especially the streamwise component. The Reynolds stresses 
near the wall gradually decrease along the airfoil. The Reynolds stresses in the wake part of the boundary layer 
gradually increase along the airfoil. The maximum Reynolds stresses occur in the inflectional part of the boundary 
layer. Notice that the boundary layer on the upper surface is evolving in an adverse pressure gradient region. In the 
DNS'*'’ and in the experiments'*'? performed over flat plates with adverse pressure gradients, the Reynolds 
stresses distributions had two peaks. One peak was near the wall region and the second peak was in the wake part of 
the boundary layer. In the current simulation, the results show only one peak in the outer part of the boundary layer 
and the second peak that existed downstream of the leading-edge separation bubble gradually decayed. There are 
two differences in flow over airfoils at high angles of attack compared to that over a flat plate in adverse pressure 
gradients. One is the existence of a laminar separation bubble near the leading edge and the other is the separation of 
the evolving boundary layer in the trailing-edge part of the airfoil. The separation bubble near the leading edge 
generates a shear layer downstream of the separation point. The flow features downstream of this region strongly 
depend on the thickness and the velocity distribution of this shear layer. Hence the status of this bubble is important 
in determining the aerodynamic properties of the airfoil. In a Large-Eddy-Simulation of flow around the 
Aerospatiale A-profile” at high angle of attack near stall, comparing with experiments, it was found that the 
treatment of the laminar separation bubble is essential in predicting the aerodynamic quantities and the 
characteristics of the boundary layer in the aft part of the airfoil. Hence, in the DNS/LES simulations the exact 
conditions - such as turbulence levels and whether the boundary layer was tripped and how it was tripped — have to 
be included. 

The contours of the average statistical quantities -- (a) kinetic energy, (b) shear stress, and (c) production — are 
plotted in Figs. 25(a-d). These quantities play important roles in RANS modeling and development. As observed in 
Fig. 24, the turbulent kinetic energy and shear stress have large values inside and immediately downstream of the 
leading-edge separation bubble. Beyond the separation bubble, the shear layer grows in thickness and the turbulent 
kinetic energy and the shear stress grow in the middle of the shear layer. Turbulent kinetic energy and shear stress 
again grow inside the wake region starting from the trailing edge. Turbulent shear stress becomes negative in the 
wake region close to the trailing edge. The turbulent production peaks near the wall and in the middle of the shear 
layer in the region immediately downstream of the separation bubble. Downstream, the production increases in the 
middle of the boundary layer and decreases near the wall. Similar to the turbulent kinetic energy and shear stress, the 
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production becomes large in the wake region behind the trailing edge. The turbulent dissipation is not resolved 
accurately with this grid sizes and hence not presented in this paper. The grid sizes in the streamwise and normal 
directions will have to be increased to resolve the turbulent dissipation accurately. 
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Figure 24. Normal and shear Reynolds stress components profiles at different stations. 
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Figure 25. Contours of the turbulent statistical quantities (a) Kinetic energy, (b) Shear stress, and (c) Production. Re, 
= 1*10°, o=15°. 
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Ill. Conclusions 


DNS were performed for flows over an NACA-0012 airfoil at a low and a moderate chord Reynolds numbers of 
Re.= 50*10°, and 1*10°, respectively, to investigate the flow dynamics, the turbulent statistical quantities, and the 
different terms in the turbulent kinetic energy balance equation. The computations were conducted at 5 degrees 
angle of attack at the low Reynolds number case and at 15 degrees for the higher Reynolds number case. The 
freestream Mach number for the low and higher Reynolds numbers cases are M = 0.4, and 0.2, respectively. The 
simulation results at the low Reynolds number showed a long laminar separation bubble in the forward part of the 
airfoil and an attached thick shear layer in the aft part of the airfoil. For this case, the computations of the drag 
coefficient revealed that the form drag contributed about 76% of the total drag and the friction drag contributed only 
24% of the drag at an angle of attack of 5 degrees. It was observed that the flow near the reattachment consists of an 
outer region that resembles a turbulent mixing layer. It was also noted that the outer shear layer that originated from 
the separation bubble persisted up to the trailing edge at this low Reynolds number. The developing wall boundary 
layer along the wall recovers slowly toward an equilibrium boundary layer. It was also observed that turbulent 
intensities and Reynolds shear stresses peak near the reattachment region and in the middle part of the shear layer. 
The budget of the terms in the turbulent kinetic energy equation showed that the roles of the turbulent diffusion, 
advection, and dissipation terms change as the boundary layer evolves over the airfoil. In the separation bubble 
region, the dissipation is small compared to the diffusion and advection terms. Dissipation becomes larger than the 
diffusion and advection terms downstream of the reattachment point 

The simulations at larger Reynolds number and at a high angle of attack showed a small laminar separation 
bubble on the order of 2% of the chord at the airfoil leading edge (x/c ~ 0.008 to 0.028). The thickness of the 
separated region is about 0.002 times the chord. The effects of the separated shear layer diminished within a short 
distance from the reattachment point. A well-developed turbulent boundary layer with a strong velocity gradient 
near the wall formed downstream of the separation bubble. As it evolves downstream, the boundary layer thickness 
increases and an inflectional profile starts to appear even before the trailing edge separation point. Eventually, the 
boundary layer separates near the trailing edge. For this case, the computations of the drag coefficient revealed that 
the form drag contributed about 92% of the total drag and the frictional drag contributed only 8% of the drag at an 
angle of attack of 15 degrees. The results show that the Reynolds stresses in the leading edge separation bubble are 
one order of magnitude higher than in the downstream part of the airfoil. The turbulent normal stresses peak in the 
maximum shear region as observed in the low Reynolds number case. There is a sharp peak in the normal Reynolds 
stresses near the wall immediately downstream of the reattachment point. Reynolds stresses exist in the outer part of 
the boundary layer. The Reynolds stresses near the wall decrease along the chord and the Reynolds stresses 
gradually increase in the outer part of the boundary layer. The grid used in the present simulation is not sufficient to 
resolve the turbulent dissipation at this Reynolds number. The grid resolutions in the streamwise and normal 
directions have to be increased in order to capture the Kolmogorov scales. 

The current DNS revealed the complex flow features that exist in the flow over an airfoil at angles of attack. The 
flow on the suction side consists of a laminar separation bubble, a turbulent reattachment, and a trailing edge 
turbulent separation. In the next step, LES simulations are planned to validate how these methods capture these 
features. The flow did not massively separate at this angle of attack of 15 degrees. Simulations should be conducted 
at high angles of attack beyond stall to obtain the flow features and the turbulent statistics at those massively 
separated cases. As discussed before, the freestream turbulence will strongly influence the existence or the size of 
the leading-edge separation bubble. Simulations are planned incorporating freestream turbulence to investigate its 
effects on the aerodynamic and the turbulent characteristics of flows over airfoil at high angles of attack. 
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